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Abstract 

The aim of the present paper is to study the effects of Hebbian learning in random 
recurrent neural networks with biological connectivity, i.e. sparse connections and 
separate populations of excitatory and inhibitory neurons. We furthermore consider 
that the neuron dynamics may occur at a (shorter) time scale than synaptic plas- 
ticity and consider the possibility of learning rules with passive forgetting. We show 
that the application of such Hebbian learning leads to drastic changes in the network 
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dynamics and structure. In particular, the learning rule contracts the norm of the 
weight matrix and yields a rapid decay of the dynamics complexity and entropy In 
other words, the network is rewired by Hebbian learning into a new synaptic struc- 
ture that emerges with learning on the basis of the correlations that progressively 
build up between neurons. We also observe that, within this emerging structure, 
the strongest synapses organize as a small-world network. The second effect of the 
decay of the weight matrix spectral radius consists in a rapid contraction of the 
spectral radius of the Jacobian matrix. This drives the system through the "edge 
of chaos" where sensitivity to the input pattern is maximal. Taken together, this 
scenario is remarkably predicted by theoretical arguments derived from dynamical 
systems and graph theory. 

Key words: Random recurrent neural networks, Hebbian Learning, Network 
structure, Chaotic dynamics 



1 Introduction 

Neural networks show amazing abilities for information storage and process- 
ing, and stimulus- dependent activity shaping. These capabilities are mainly 
conditioned by the mutual coupling relationships between network structure 
and neuron dynamics. Actually, learning in neural networks implies that ac- 
tivity guides the way synapses evolve; but the resulting connectivity structure 
in turn can raise new dynamical regimes. This interaction becomes even more 
complex if the considered basic architecture is not feed-forward but includes 
recurrent synaptic links, like in cortical structures. Understanding this mutual 
coupling between dynamics and topology and its effects on the computations 
* Corresponding author: hugues.berry@inria.fr 
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made by the network is a key problem in computational neuroscience, that 
could benefit from new approaches. 

In the related field of dynamical systems interacting via complex coupling 
networks, a large amount of work has recently focused on the influence of 
network topology on global dynamics (for a review, see [9]). In particular, 
much effort has been devoted to understanding the relationships between node 
synchronization and the classical statistical quantifiers of complex networks 
(degree distribution, average clustering index, mean shortest path, modular- 
ity...) [21,35,32]. The main idea was that the impact of network topology on 
the global dynamics might be prominent, so that these structural statistics 
may be good indicators of the global dynamics. This assumption proved how- 
ever largely wrong so that some of the related studies yielded contradictory 
results [35,27]. Actually, synchronization properties cannot be systematically 
deduced from topology statistics but may be inferred from the spectrum of the 
network [3]. Accordingly, many studies have considered diffusive coupling of 
the nodes [22]. In this case, the adjacency matrix has real nonnegative eigen- 
values, and global properties, such as stability of the synchronized states [4], 
can easily be inferred from its spectral properties. 

In this perspective, neural networks can be considered as mere examples of 
these complex systems, with the particularity that the dynamics of the net- 
work nodes (neurons) depends on the network links (synaptic weights), that 
themselves vary over time as a function of the node dynamics. Unfortunately, 
the coupling between neurons (synaptic weight) is rarely diffusive, so that the 
corresponding matrix is not symmetric and may contain positive and negative 
elements. Hence the mutual coupling between neuron dynamics and network 
structure remains largely to be understood. 

Our general objective is to shed light on these interactions in the specific 
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case of random recurrent neural networks (RRNNs). These network models 
display a rich variety of dynamical behaviors, including fixed points, limit 
cycle oscillations, quasiperiodicity and deterministic chaos [16], that are sus- 
pected to be similar to activity patterns observed in the olfactory bulb [41,18]. 
It is known that the application of biologically-plausible local learning rules 
(Hebbian rules) reduces the dynamics of chaotic RRNNs to simpler attractors 
that are specific of the learned input pattern [13]. This phenomenon endows 
RRNNs with associative memory properties, but remains poorly understood. 
Our previous work showed that the evolution of the network structure during 
learning can be tracked in numerical simulations via the classical topological 
statistics from "complex networks approaches" [6]. More recently, we devised 
a mathematical framework for the effects of Hebbian learning on the dynam- 
ics, topology and some functional aspects of RRNNs [40]. This theoretical 
approach was shown to explain the effect of learning in a "canonical" RRNN, 
i.e. a completely connected network where a neuron projects both excitatory 
and inhibitory synapses. However, this network type remains rather poorly 
realistic from a biological point of view. 

The aim of the present paper is thus to study the effects of a more biological 
connectivity. In particular, we segregate the neurons into two distinct popula- 
tions, namely excitatory (projecting only excitatory synapses) and inhibitory 
(projecting only inhibitory synapses) neurons. Furthermore, the network is 
sparsely connected and the overall connectivity parameters are fixed to em- 
ulate local circuitry in the cortex. We show that the application of Hebbian 
learning leads to drastic changes in the network dynamics and structure. We 
also demonstrate that the mathematical arguments mentioned above remain 
a very useful unifying framework to understand the effects of learning in this 
system. 
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2 The model 



2. 1 Connectivity 

We consider networks with a total oi N = 500 neurons and random connectiv- 
ity. Each neuron is either inhibitory (with probability pi) or excitatory (with 
probability Pe = 1 — Pi) and projects to p c N randomly chosen postsynaptic 
neurons (independently of their excitatory or inhibitory nature). Probabilities 
are taken uniform on [0, 1] for the network connectivity p c and fraction of 
inhibitory neurons pi. In the present study, we fixed pi and p c so as to ac- 
count for the neural circuitry of a typical neocortical column. Hence, we used 
pi = 0.25 [30] and p c = 0.15 [33,29]. 

The initial weight of each synapse between a postsynaptic neuron i and a 
presynaptic neuron j, is drawn at random, according to a Gamma dis- 

tribution, whose parameters depend on the nature of the presynaptic neuron 
j. If j is inhibitory, wjp ~ Gamma(— p, w /rii, cr w /ni), where Gamma(m, s) de- 
notes the Gamma distribution with mean m and standard deviation s, and 
n-i = PiPcN. If j is excitatory, then wjp ~ Gamma(/x ti ,/n e , a w /n e ) where 
n e = PePcN . Using Gamma distributions (instead of Gaussian ones, for in- 
stance) allows to ensure that inhibitory (excitatory) neurons project only neg- 
ative (positive) synapses, whatever the values of ji w and a w . Thanks to the 
normalization terms (n e and rij), the total excitation received by a postsynap- 
tic neuron is on average equal to the total inhibition it receives. Hence, in their 
initial setups (i.e. before learning) our networks are guarantied to conserve the 



5 



excitation/inhibition balance (on average). 



2.2 Dynamics 

We consider firing-rate neurons with discrete-time dynamics and take into ac- 
count that learning may occur on a different (slower) time scale than neuron 
dynamics. Synaptic weights are thus kept constant for r > 1 consecutive dy- 
namics steps, which defines a "learning epoch" . The weights are then updated 
and a new learning epoch begins. We denote by t > the update index of 
neuron states (neuron dynamics) inside a learning epoch, while T > 1 indi- 
cates the update index of synaptic weights (learning dynamics). 
Let x\ T \t) G [0, 1] be the mean firing rate of neuron i, at time t within the 
learning epoch T. Let be the matrix of synaptic weights at the T-th 

learning epoch and £ the vector ((,i)f =1 . Then the discrete time neuron dy- 
namics (1) writes: 



Here, / is a sigmoidal transfer function (f(x) = 1/2 (1 + tanh(prr))). The 
output gain g tunes the nonlinearity of the function and mimics the excitability 
of the neuron. is a (time constant) external input applied to neuron i and 

(T) 

the vector £ is the "pattern" to be learned. W-j represents the weight of the 
synapse from presynaptic neuron j to postsynaptic neuron i during learning 
epoch T. Finally, at the end of one learning epoch, the neuron dynamics indices 





are reset: x\ 



(T+l) 



(0) = x\ 1 >(r),\/ i . 
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2. 3 Learning 



In the present work, we used the following Hebbian learning rule (see [40] for 
justifications): 

= AWf > + s^mPmPH (mf >) . (2) 

where a is the learning rate, Sj — +1 if j is excitatory and —1 if it is inhibitory 
and H denotes the Heaviside step function (H(x) — if x < 0, 1 else). The 
first term in the right-hand side (RHS) member accounts for passive forgetting, 
i.e. A G [0, 1] is the forgetting rate. If A < 1 and or nij = (i.e. the pre- 
or postsynaptic neurons are silent, see below), eq.(2) leads to an exponential 
decay of the synaptic weights (hence passive forgetting). Another important 
consequence of this rule choice is that if A < 1, the weights are expected to 
converge to stationary values. Hence A < 1 also allows avoiding divergence of 
the synaptic weights. Note that there is no forgetting when A = 1. 
The second term in the RHS member of eq. (2) generically accounts for activity- 
dependent plasticity, i.e. the effects of the pre- and postsynaptic neuron firing 
rates. In our model, this term depends on the history of activities through the 
time-average of the firing-rate: 

mP = -±xP(t)-d t , (3) 
T t=i 

where di G [0, 1] is a threshold that we set to di = 0.10, \/i in the present study. 
A neuron i will thus be considered active during learning epoch T whenever 
mP > (i.e. whenever its average firing rate has been > 10% of the maximal 
value), and silent else. 

Note that definition (3) actually encompasses several cases. If r = 1, weight 
changes depend only on the instantaneous firing rates, while if r ^> 1, weight 
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changes depend on the mean value of the firing rate, averaged over a time 
window of duration r in the learning epoch. In many aspects the former case 
can be considered as genuine plasticity, while the latter may be related to 
meta-plasticity [1]. In this paper, we used r = 10 4 . Finally, weights cannot 
change their sign. Note however that this setup does not have a significant 
impact on the present results. 

3 Results 

3.1 Spontaneous dynamics 

We first present simulation results on the spontaneous dynamics of the sys- 
tem, i.e. the dynamics eq.(l) in the absence of learning. The phase diagrams 
in fig. 1 locate regions in the parameter space for which chaotic dynamics are 
observed in simulations. Each panel shows the isocurve L 1 = (where Li is 
the largest Lyapunov exponent) that represents the boundary between chaotic 
(Li > 0) and non chaotic (Iq < 0) dynamics. 

It is clear from this figure that chaotic dynamics are found for large parts of 
the parameter space. Generally speaking, chaotic behaviors disappear when 
the average weight fj, w increases, which may be related to an increase of the 
average neuron saturation. A more surprising observation is that chaotic dy- 
namics tends to disappear when the gain of the transfer function g is increased. 
This behavior is in opposition to the behavior observed with classical random 
recurrent networks with homogeneous neurons (where each neuron has both 
excitatory and inhibitory projections). In the latter models (and even in re- 
lated two-populations models, see [14]), chaotic dynamics usually appear for 
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increasing values of g (see e.g. [11]). 

This is a very interesting property of the spontaneous dynamics in our model, 
whose understanding is however out of the scope of the present paper and is 
left for future work. In the framework of the present study, these phase dia- 
grams mainly allow locating suitable parameters for the initial conditions of 
our networks. We wish the initial dynamics to provide a wide range of possi- 
ble dynamical regimes, a large (KS) entropy and self-sustaining dynamics. For 
these reasons, we set our initial dynamics inside the chaotic region, and fix 
/j, w = 50, a w = 1.0, g — 10 and N = 500. The initial weight distribution will 
thus consist in a Gamma distribution with effective average —2.67 and s.d. 
0.053 for inhibitory synapses, and 0.89 and 0.018, respectively, for excitatory 
ones. 

3.2 Structure modifications 

In this section, we study the changes in the network structure induced by the 
learning rule described by eq.(2). 

3.2.1 Adjacency matrix 

The adjacency matrix of a graph gives information about its connectivity and 
can be extracted from the weight matrix W by thresholding and binarization. 
We applied a simple relative thresholding method that consists in keeping only 
the absolute values of the 9 percent highest weights (again, in absolute value) 
from the nonzero connections in W. Hence gradual decrease of 6 enables to 
progressively isolate the adjacency network formed by the most active weights 
only. The resulting matrix is then binarized and symmetrized, yielding the 
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adjacency matrix A{6) whose elements a,ij(0) indicate whether i and j are 
connected by a synapse with a large (> 9) weight (either inhibitory or exci- 
tatory), compared to the rest of the network. We limit the range of 6 values 
to ensure that not more that 10% of the neurons get disconnected from the 
network by the thresholding process. 

To characterize the topology of these matrices, we computed the two main 
quantifiers used in "complex networks" approaches, namely the clustering in- 
dex and the mean shortest path (see [40] for formal definitions). The clustering 
index C is a statistical quantifier of the network structure and reflects the de- 
gree of "cliquishness" or local clustering in the network [47]. It expresses the 
probability that two nodes connected to a third one are also connected to- 
gether and thus can be interpreted as the density of triangular subgraphs in 
the network. The mean shortest path (MSP) is the average, over all noniden- 
tical neurons pairs of the smallest number of synapses one must cross 

to reach % from j. 

These two quantifiers are usually informative only when compared to similar 
measures obtained from reference random networks [47]. Here, to build ref- 
erence networks, we start with the weight matrix at learning epoch T, 
and rewire it at random but preserving the inhibitory/excitatory nature of 
the neurons. Hence for each element W^J\ we choose (uniformly) at random 

(T) 

another element with the same sign, and exchange their values. We then 
compute the clustering index and mean shortest path of the resulting rewired 
network, and average the obtained values over 15 realizations of the rewiring 
process, yielding the reference values C ran d{0) and M S P ran d{9) . 
Figure 2A & B show simulation results for the evolution of C^ T \6) / C^ nd (9) 
and MSP (T \9)/MSP^ d (0) during learning. The distribution of the initial 
weights over the network being random, the resulting adjacency matrix (9) 
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is essentially a random network, i.e. one expects {9) / C^ nd {9) 1 and 
MSP^(9)/MSPiH d (9) « 1, V0. This is confirmed in fig. 2: during the ap- 
proximately first 100 learning epochs, both network statistics remain around 
1. 

The situation however changes for longer learning epochs. For T > 100, the 
relative MSP remains essentially 1 for all thresholds 9 (less than 4% varia- 
tion, fig. 2B). Hence, the average minimal number of synapses linking any two 
neurons in the network remains low, even when only large synapses are con- 
sidered. Conversely, the clustering index (fig. 2 A) increases at T > 100 for the 
stronger synapses and reaches a stable value that is up to almost two twofold 
the value found in the reference random networks. Hence, if one considers the 
strong synapses at long learning epochs, the probability that the neighbors of 
a given neuron are themselves interconnected is almost twofold higher than if 
these strong synapses were laid at random over the network. In other terms, 
the learning rule yields correlations among the largest synapses at long learn- 
ing epochs. 

In the literature related to "complex networks", networks with a larger cluster- 
ing index but a similar MSP with respect to a comparable reference random 
network, are referred to as small-world networks. Hence, the learning rule 
eq.(2) organizes strong synapses as a "small- world" network. 
Emerging experimental evidence shows that numerous brain anatomical and 
functional connectivity networks at several length scales indeed display a com- 
mon small- world connectivity (for a recent review, see [5]). These include 
quantifications of the physical [37] or functional [8] connectivity of neuronal 
networks grown in vitro; quantifications of the anatomical connectivity of 
Caenorhabditis elegans full neural system [47] or, at larger scale, cortico- 
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cortical area connectivity maps in macaque, cat [43] and more recently hu- 
man [23]; and quantitative studies of functional human brain connectivity 
based on MEG [44], EEG [34] or fMRI data [2,17]. Current hypothesis for the 
frequent observation of small-world connectivity in real biological networks 
state that they may result from natural evolution (small-world networks are 
generally thought to minimize wiring length while preserving low energy costs 
and high local integration of the neurons [5,28]). 

An alternative hypothesis could be that small-world networks are emerging 
properties of neural networks subject to Hebbian learning. In favor of this 
possibility, small-world connectivity has recently been shown to arise sponta- 
neously from spiking neuron networks with STDP plasticity and total connec- 
tivity [38] or with correlation-based rewiring [31]. Hence our present findings 
tend to strengthen this hypothesis. 

However, these kinds of interpretation should be taken with great care. For 
instance, it is easy to find learning rules similar to eq.(2) or areas in the pa- 
rameter space, for which the network, even at long learning times, only slightly 
deviates from its initial random organization (see e.g. [40]). Hence, emergence 
of small-world connectivity, even in computational models of neural networks 
(i.e. not to speak about real neural networks), may be restricted to certain 
areas of the learning rule parameter space. 

More importantly, these indicators in fact give no obvious clue about the mu- 
tual coupling between global dynamics and the network structure. Hence, in 
our case at least, the classical statistics of the "complex networks" do not pro- 
vide causal explanation for the dynamical effects of learning. For instance, it 
does not help understand why dynamics complexity systematically decreases 
during learning. However, the adjacency matrix is not the only viewpoint from 
which the network structure can be observed (see [40] for a discussion). In the 
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following, we examine what information can be obtained if the structure is 
observed at the level of the Jacobian matrices. 

3.2.2 Jacobian matrices. 

Denote by F the function F : -> R N such that F;(x) = /(a*). In our case, 
the components of the Jacobian matrix of F at x, denoted by -DF X are given 
by: 

dF N 

= f'(E W ^k + tiWij = f(u t )W ir (4) 
ax i fc=i 

Thus it displays the following specific structure: 

DF X = A(u)W, (5) 

with: 

Atf(u) = f'i^Sij. (6) 

Note that DF X depends on x, contrarily to W. Generally speaking, -DF X gives 
the effects of perturbations at the linear order. To each Jacobian matrix -DF X 
one can associate a graph, called "the graph of linear influences". To build 
this graph, one draws an oriented link j — > i iff d/K) ^ q. The link is 
positive if > and negative if < 0. A detailed presentation of the 
properties of the graph of linear influences can be found in [11,40]. We just 
recall here that this graph contains circuits or feedback loops. If e is an edge, 
we denote by o(e) the origin of the edge and t(e) its end. Then a feedback 
loop is a sequence of edges e±, ...,e£ such that o(e i+1 ) = t(ei), V? = l...k — 1, 
and t(ek) = o(ei). A feedback loop is said positive (negative) if the product of 
its edges is positive (negative). 

In general, positive feedback loops are expected to promote fixed-point stabil- 
ity [25] whereas negative loops usually generate oscillations [45,20]. In a model 
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such as eq.(l) the weight of a loop k±, k 2 , . . . k n , k\ is given by n™=i ^h +1 kif'(uki), 
where k n+ i = k±. Therefore, the weight of a loop is the product of a "topolog- 
ical" contribution (IJJLi Wfc i+ ifcj) and a dynamical one (njLi f'( u ki))- 
We measured the evolution of feedback loops during learning via the weighted- 
fraction of positive circuits in the Jacobian matrix, that we defined as 

where (resp. <r~^) is the sum of the weights of every positive (resp. neg- 

ative) feedback loops of length n in the Jacobian network at learning epoch T. 
Hence G [0, 1]. If its value is > 0.5, the positive feedback loops of length 
n are stronger (in total weight) than the negative ones in the network. We 
computed the weighted-fraction of positive feedback loops for length n = 2 
and n = 3 (i.e. and R^). 

The evolutions of R 2 and R\ are presented in fig. 3A. During the first ~ 20 
learning epochs, the time course of these quantities are highly noisy (and the 
corresponding standard deviation very large), so that we could not interpret 
them conclusively. However, a T 25 learning epochs, R^ stabilizes to val- 
ues < 0.5 (i?^ ~ 0.47), indicating a slight imbalance in favor of negative 
feedback loops over positive ones. According to the above theoretical con- 
siderations, this indicates a trend toward complex oscillatory dynamics. This 
viewpoint may be considered another perspective to explain the initial chaotic 
dynamics. Note however that the initial imbalance in circuits of length-3 is 
much more modest, R^ ss 0.497. 

(T) 

When 25 < T < 50, R 2 increases and converges to w 0.50. A dynamical 
interpretation would be that the corresponding dynamics attractors become 
progressively less chaotic and more periodic. This is exactly the behavior ob- 
served in the simulations (see fig. 5B). Hence in spite of the huge fluctuations 
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observed at the beginning, the study of the feedback loops in the Jacobian 
matrix offers a useful interpretation to the reduction of dynamics induced by 
learning at short learning epochs. 

Upon further learning, and remain constant at 0.5 up to T pa 100 
learning epochs. Thus, these quantities do not detect variations in the balance 
between positive and negative feedback loops for 50 < T < 100. However, 
at longer times (T > 100), R$p and R^ both increase abruptly and rapidly 

(T) (T) 

reach pa 0.62 for R 2 and ~ 0.56 for R 3 . Hence, at long learning epochs, the 
system switches to a state where positive feedback loops hold a significantly 
larger weight as compared to negative ones. Note that the time course of these 
indicators for T > 25 closely follows the time course of the relative clustering 
index (fig. 2A). The causal relation between these two phenomena is however 
not obvious. 

Because of the particular form of the Jacobian matrix in our system, the sign 
of a feedback loop is given by the sign of the weights along it (see above). We 
thus proceeded (fig. 3B) to the computation of the evolution of the weighted- 
fraction for feedback loops computed in W, i.e. we compute here the weight of 
a feedback loop e±, . . . , as the product of the synaptic weights of its edges, 
thus independently of the neuron state. The evolution of the weighted-fraction 
of positive feedback loops in W does not account for the initial imbalance ob- 
served in the feedback loops of DF. However, its evolution at long times is 
remarkably identical to that measured in DF. Thus, the weighted-fraction of 
positive feedback loops in W is able to account for at least part of the evo- 
lution of the dynamics and represents a link between purely structural and 
purely dynamical aspects. However, more information can be extracted by a 
more dynamical approach. 



15 



3.3 Dynamical perspective 



Starting from spontaneous chaotic dynamics, application of the Hebbian learn- 
ing rule (2) in our sparse two-populations model leads to dynamics simplifica- 
tion, as in the case of completely-connected, one-population random recurrent 
neural networks [13]. Figure 5B shows the network- averaged neuron dynamics 
obtained at different learning epochs. The dynamics, initially chaotic (T = 1), 
gradually settles onto a periodic limit cycle (T = 270), then on a fixed point 
attractor at longer learning epochs (see e.g. T = 290 in this figure). This 
evolution of the global dynamics is a typical example of the reduction of the 
attractor complexity due to the mutual coupling between weight evolution and 
neuron dynamics. 

In [40], we developed a theoretical approach derived from dynamical systems 
and graph theory and evidenced that it explains this reduction of complexity 
in homogenous (single population) recurrent neural networks. We shall show 
thereafter that it also provides a useful framework for the present model. Be- 
low, we first summarize the main results obtained from this mathematical 
analysis (for details, see [40]). 

3.3.1 Main theoretical results 

The first prediction of our approach is that Hebbian learning rules contract 
the norm of the weight matrix W. Indeed, we could compute the following 
upper bound: 

IIW^HA^IIW^II + ^^C, (8) 

where |||| is the operator norm (induced e.g. by Euclidean norm) and C a 
constant depending on the details of the rule. Hence the major effect of the 
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learning rule is expected to be an exponentially fast contraction of the norm 
(or equivalently the spectral radius) of the weight matrix, which is due to the 
term A, i.e. to passive forgetting (A < 1). 

The next prediction concerns the spectral radius of the Jacobian matrix. In- 
deed, one can derive a bound for the spectral radius of DF^: 

\^( X )\<m^f\vP)\\W^\\. (9) 

This equation predicts a contraction of the spectrum of DF^p that can arise 
via two effects: either the contraction of the spectrum of and/or the 

decay of maxj f'(ui), which arises from saturation in neuron activity. Indeed, 
f'iui) is small when Xi is saturated to or 1, but large whenever its synaptic 
inputs are intermediate, i.e. fall into the central part of the sigmoid f{ui). We 
emphasize that when A = 1, and u^ T ) can diverge and lead maxj f'{v!*p) 
to vanish. Hence the spectral radius of the Jacobian matrix can decrease even 
in the absence of passive forgetting. In all cases, if the initial value of |/4 T ^(x)| 
is larger than 1, eq.(9) predicts that the spectral radius may decrease down to a 
value < 1. Note that in discrete time dynamical systems the value \fj\ (x)| = 1 
locates a bifurcation of the dynamical system. 

(T) 

According to our present setting, the largest Lyapunov exponent, L\ depends 
on the learning epoch T. We were able to show that: 

4 T) < log(||W^>||) + (log(max/'(^))) , (10) 

where (log(max; f{uj)))^ denotes the time average of log(maxj /'(u,)), in 
the learning epoch T (see [40] for formal definitions). The second term in the 
RHS member is related to the saturation of neurons. The first one states that 
will decrease if the norm of the weight matrix ||W^ T ^|| decreases during 
learning, resulting in a possible transition from chaotic to simpler attractors. 
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Let u\ \t) = Y,f=i Xj '(t) + the local field (or membrane potential) 
of neuron i at dynamics step t within learning epoch T. Our theoretical work 
also showed that provided A < 1, the vector u = (ui)f =1 converges to a fixed 
point as T — > +00: 



Therefore, the asymptotic local field is predicted to be the sum of the input 
pattern plus an additional term H*- 00 -*, which accounts for the history of the 
system and can be weak or not, depending on the exact learning rule and 
system history. 

In [40] , we studied the effects of Hebbian learning in a completely connected 
(p c = 1) one-population network (i.e. where each neuron can project inhibitory 
(negative) and excitatory (positive) synapses) and showed that these analyt- 
ical arguments explain and describe results of the system simulation with a 
very good accuracy. 

While the model studied in the present work is much more compatible with 
our knowledge of biological neural networks, it is very different from the model 
studied in [40]. In the present model, the connectivity is (severely) sparse and 
the neurons are segregated in two distinct groups, with distinct synaptic prop- 
erties. Furthermore, the learning rule eq.(2) is also more complex. Hence, it 
is not clear whether the above theoretical arguments account for the current 
case. In particular, these arguments mainly provide upper bounds, whose qual- 
ity is not guarantied. In the following sections, we present simulation results 
about the influence of learning on the network dynamics and function, using 
our theoretical framework as an oracle. 




(00) 



| + H (oo) , 



(11) 



where: 




(12) 
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3.3.2 Dynamics evolution in the sparse 2-populations model. 

Figure 4 shows the evolution of the spectral radius of W, |sf^| for A = 0.90 or 
0.99 in simulations of our sparse two-populations model with dynamics eq.(l) 
and learning rule eq.(2). Let SjP be the eigenvalues of W^ T \ ordered such that 

l s i T ' ) | > l^l > • • • > 4^ > Since Is^j, the spectral radius of W^ T \ is 

smaller than ||W^|| one has from eq.(8): 

I^I^A-HWWll + ^^a (13) 

It is clear from this figure that in both cases the spectral radius decreases ex- 
ponentially fast, with a rate that is very close to the prediction of the theory 
(i.e. oc \ T ). Hence, the decay predicted by our analytical approach (eq.8) is 
obviously observed in the simulations. Note that the clear trend in the simu- 
lation results for a decay proportional to A T , even tells us that the bound in 
(13) is indeed very good. 

Figure 7 shows (among other curves) the evolution of |/4 T ^(x)| (dashed thin 
line). This figure confirms that the theoretical prediction about the decay of 
|/4 T ^(x)| (eq.9) is also valid for this model. Hence, eq.(9) opens up the pos- 
sibility that learning drives the system through bifurcations. This aspect is 
studied below (section 3.3.3). 

We now turn to directly study how the attractor complexity changes dur- 
ing learning. This information is provided by the computation of the largest 
Lyapunov exponent. Note that another canonical measure of dynamical com- 
plexity is the Kolmogorov-Sinai (KS) entropy which is bounded from above by 
the sum of positive Lyapunov exponents. Therefore, if the largest Lyapunov 
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exponent decreases, the KS entropy decreases as well. 

(T) 

Figure 5A shows the evolution of L\ during numerical simulations with dif- 
ferent values of the passive forgetting rate A. Its initial value (L^ pa 0.94) is 
positive (we start our simulations with chaotic networks). As predicted by our 
theoretical approach (eq.10), the Hebbian learning rule eq.(2) leads to a rapid 
decay of The decay rate is indeed close to log(|| ||) for intermediate 
learning epochs, in agreement with the upper bound of eq.(10). Hence 
quickly shifts to negative values, confirming the decrease of the dynamical 
complexity that could be inferred from visual inspection of fig. 5B. 
One can also consider individual neuron activities. Figure 6 shows the evolu- 
tion of the local field u during learning. Clearly, the initial values are random, 
but the local field (thin gray line) shows a marked tendency to converge to 
the input pattern (thick dashed line) after as soon as 60 learning epochs. At 
T = 180, the convergence is almost complete. Hence this behavior once again 
conforms to the theoretical predictions eq.(ll), with £ 3> H*- 00 - 1 . In the results 
presented in this figure, we pursue the simulation up to T = 200, at which 
point we remove the pattern from the network, i.e. we set & = 0, Mi (fig. 2D). 
As a result, u looses its alignment from the pattern and presents a noisy as- 
pect (note that each vector in the figure has been normalized to [0, 1]). This 
behavior is once again in agreement with the theoretical predictions of eq.(ll), 
which indicates that (u)^ = H^°°^ upon pattern removal. 
To conclude, we have shown here that Hebbian learning in our system leads 
to a decrease of the attractor complexity and entropy that can be induced by 
passive forgetting and/or an increased level of saturation of the neurons. This 
corresponds in details to the scenario predicted by our mathematical analysis. 
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3.3.3 Functional consequences 



The former sections dealt with the effects of Hebbian learning on the structure 
and dynamics of the network. We now turn to its influence on the network 
function. The basic function of RRNNs is to learn a specific pattern £. In 
this framework, a pattern is learned when the complex (or chaotic) dynamics 
of the network settles onto a periodic oscillatory regime (a limit cycle) that 
is specific of the input pattern. This behavior emulates putative mechanisms 
of odor learning in rabbits that have been put forward by physiologists such 
as W. Freeman [18,19]. An important functional aspect is that removal of 
the learned pattern after learning should lead to a significative change in the 
network dynamics. We now proceed to an analysis of this latter property. 
The removal of £ is expected to change the attractor structure and the average 
value of any observable (f> (though with variable amplitude). Call A^ T ^ [(f)] the 
variation of the (time) average value of induced by pattern removal. If the 
system is away from a bifurcation point, removal will result in a variation of 
A^ T ) [(f)] that remains proportional to 

On the opposite, close to a bifurcation point this variation is typically not 
proportional to £ and may lead to drastic changes in the dynamics. From 
the analysis above, we therefore expect pattern removal to have a maximal 
effect at "the edge of chaos" , namely when the (average) value of the spectral 
radius of DF X is close to 1. Call and v fc the eigenvalues and eigenvectors 
of W (T) A(u*( T )), ordered such that |A^| < |Ajv_i| < |Ai| < 1. In the case 
where the dynamics has converged to a stable fixed point u*( T ) (namely, when 
L\ < 0, see e.g. fig. 5), our theoretical work predicted that: 

a( t) [u] = _ £ <p^|) Vfe (14) 

k=l 1 - A k 
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where ( , ) denotes the inner product. As a matter of fact, the RHS term 
diverges if Ai = 1 and if (vi,£) 7^ 0. From this analysis, we therefore expect 
pattern removal to have a maximal effect at "the edge of chaos" , namely when 
the value of the spectral radius of DF X is close to I. 

To study the effects of pattern removal in our model, we monitored the quan- 
tity: 



1 N 9 
A (T) [A] = ^E((A,(u)) (r) -(A,(u')) (T) ) (15) 

\ i=l 

that measures how neuron excitability is modified when the pattern is re- 
moved. The evolution of A( T )[A] during learning with rule eq.(2) is shown on 
fig. 7 (thick full lines) for two values of the passive forgetting rate A. A^ T ^ [A] 
is found to increase to a plateau, and vanishes afterwards. Interestingly, com- 
parison with the decay of the leading eigenvalue of the Jacobian matrix, ^ 
(thin full lines) shows that the maximal values of A( T )[A] are obtained when 

is close to 1 and the largest Lyapunov exponent L\ close to 0. 
Hence, these numerical simulations are in agreement with the theoretical pre- 
dictions that Hebbian learning drives the global dynamics through a bifurcation, 
in the neighborhood of which sensitivity to the input pattern is maximal. Note 
that this property is obtained at the frontier where the chaotic strange attrac- 
tor begins to destabilize (|//i| = 1), hence at the so-called "edge of chaos". This 
particular dynamical regime, at the frontier between order (periodical or fixed 
point regimes) and disorder (chaos), has already be reported to be particularly 
suitable for recurrent neural networks, especially when computational power 
is considered [42,7]. The present results show that it is the optimal regime for 
the sensibility to the input pattern in our model. Whether this also implies 
improved or optimal computational performance remains however to be tested 
and will be the subject of future works. 
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It must finally be noticed that our theory predicts that pattern sensitivity 
should be maximal when is close to one. But several aspects of our simu- 
lation results are not accounted for by this theory. For instance, fig. 7 shows 
that \/jLi\ approaches 1 at several learning epochs. This is related to the "Arnold 
tongue" structure of the route to chaos. However, pattern sensibility is max- 
imal only for the last episode, and almost zero for the former ones. This 
behavior is still unclear and will be the subject of future works. 



4 Conclusion and future works 

To conclude, we have shown in this work that Hebbian learning eq.(2) has 
important effects on the dynamics and structure of a sparse two-populations 
RRNN. The forgetting part of the learning rule contracts the norm of the 
weight matrix. This effect, together with an increase in the average saturation 
level of the neurons, yields a rapid decay of the dynamics complexity and en- 
tropy. In other words, the network forgets its initial synaptic structure and is 
rewired by Hebbian learning into a new synaptic structure that emerges with 
learning and that depends on the whole history of the neuron dynamics. We 
have shown that the strongest synapses organize within this emerging struc- 
ture as a small-world connectivity. The second effect of the decrease of the 
weight matrix and of the increased neuron saturation consists in a rapid con- 
traction of the spectral radius of the Jacobian matrix. This leads the system 
to the edge of chaos, where sensitivity to the input pattern is maximal. This 
scenario is remarkably predicted by the theoretical arguments we developed 
in [40]. 
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In the presented simulations, most of the effects are mediated by the passive 
forgetting term. We believe that this term is not unrealistic from a biological 
point of view. Indeed, synaptic plasticity at the single synapse level is not per- 
manent and some studied reported durations of 20 mn [46] or even 20 sec [10]. 
This would be accounted for in our model by A <C 1. 

Nevertheless, most studies about long-term plasticity have evidenced longer 
cellular memory time constants, ranging from hours to days [24,36,15], which 
would correspond in our model to higher A values. Note however that according 
to our mathematical analysis, most of the effects reported here are expected to 
occur even without passive forgetting (i.e. with A = 1), provided the learning 
rule increases the average saturation of the neurons. In previous studies, we 
have considered such Hebbian learning rules devoid of passive forgetting but 
provoking increasing average saturation levels of the neurons. Numerical sim- 
ulations have clearly evidenced a reduction of the attractor complexity during 
learning, in agreement with this suggestion [6,39]. 

Future works will focus on the study of more detailed biological learning 
rules (heterosynaptic LTD, synaptic rescaling). We will also consider activity- 
dependent synaptic turnover (pruning/sprouting). Indeed albeit an overlooked 
phenomena for several decades, synaptic (or at least dendritic) turnover is now 
recognized as an important part of cortical networks, even in the adult (see 
e.g. [26]). Finally, one important problem with the application of RRNNs as 
artificial neural networks, is that it is very difficult to determine when to stop 
the learning process. Our results show that the effect of an input pattern is 
maximal at those learning epochs when the system is close to a bifurcation, 
but much more modest for shorter and longer learning times. One interesting 
development would thus consist in trying to find learning rules or settings 
that would guaranty that the system remains close to the edge of chaos, even 
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at long learning times. As an attractive possibility, the plasticity of intrinsic 
properties [12] could allow the network to stabilize its activity in this region. 
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Figure 1. Phase diagram for the spontaneous dynamics eq.(l). The full line rep- 
resents the boundary between chaotic and non chaotic dynamics (i.e. the isocurve 
L\ = where L\ is the largest Lyapunov exponent). Shown are projection in (A) 
the (g,/j, w ) plan with a w = 1.0 or (B) the (a w ,fi w ) parameter plan with g = 10.0. 
Other parameters were: & = 0.0 Vi = 1 . . . N, p c = 0.15, pi = 0.25 and ./V = 500. 
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Figure 2. Evolution of the normalized structural statistics during learning with 
rule eq.(2). Values are averages over 20 different realizations of the network (ran- 
dom initial firing rates and synaptic weights). The values of the threshold 9 are, 
from bottom to top in each panel, 100%, 87%, 73%, 60% and 47%. (A) Nor- 
malized clustering index C^ T \6)/C^^ d (9). (B) Normalized mean-shortest path 
M S p( T \6) I M S Pj-and^) • Bars are ±1 standard deviation. Other parameters were: 
A = 0.90, a = 5 x 10~ 3 , g = 10, & = O.OlOsin (2m/ N) cos (8m /N) Vi = 1 . . . N, 
fx w = 50, a w = 1.0 and N = 500. 
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Figure 3. Evolution of the weighted-fraction of positive feedback loops Rn for loops 
in DF (A) and W (B) and circuit length n = 2 (thick line) and n = 3 (dotted line). 
Values are averages over 20 different networks using A = 0.90. Standard deviations 
are omitted for readability purpose. See text for definitions. All other parameters 
as in fig. 2. 
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Figure 4. Contraction of the spectral radius of W. The evolution during learning of 
the norm of W largest eigenvalue, \s^\ (thick lines) is plotted on a log-lin scale for 
A = 0.90 (bottom) or 0.99 (top). Each curve is an average over 20 realizations with 
different initial conditions (initial weights and neuron states). For clarity, standard 
deviations are omitted but are always < 20% of the average. Dashed thin lines are 
plots of exponential decreases with equation g(T) oc X T . All other parameters as in 
fig- 2. 
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Figure 5. Reduction of the dynamics complexity from chaotic to periodic and 
fixed point. (A) Evolution of the largest Lyapunov exponent L\ (full thick lines) 
for A = 0.90 (bottom) or 0.99 (top). Each value is an average over 20 realiza- 
tions with different initial conditions (initial weights and neuron states). The thin 
dashed lines illustrate decays as g(T) oc Tlog(A). (B) Examples of network dy- 
namics when learning is stopped at (from bottom to top) T = 1 (initial condi- 
tions), 270 or 290 and for A = 0.99. These curves show the network-averaged state 
(x( T \t)') = 1/N YliLi x T\t) an d are shifted along the y-axis for readability. The 
height of the vertical black bar represents an amplitude of 0.1. All other parameters 
are as in fig. 2. 
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Figure 6. Alignment of the local field u = Wx+£ (thin gray line) with the input pat- 
tern £ (thick dashed line). Snapshot are presented at T = 1 (A, initial conditions), 
T = 60 (B), T = 180 (C) and T = 200 with pattern removed (D) learning epochs. 
Each curve plots averages over 20 realizations (standard deviations are omitted for 
comparison purposes), and every vector has been normalized to [0,1] for clarity. 
A = 0.90 and all other parameters as in fig. 2 
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Figure 7. Network sensitivity to the input pattern is maximal close to a bifurca- 

(T) 

tion. The evolution of the average value of the spectral radius of -DFx (thin full 
line) is plotted together with the sensitivity measure A( T )[A] (thick full line) for 
A = 0.90 (A) or 0.99 (B). The panels also display the corresponding evolution of 
the largest Lyapunov exponent L\, plotted as 1.0 + L\ for obvious comparison pur- 
pose (thin dashed line) . The horizontal dashed-dotted lines locates y = 1 . The values 
of A( T )[A] are normalized to the [0, 1] range for comparison purposes. Each value 
is an average over 20 realizations (standard deviations are omitted for clarity). All 
other parameters were as in fig. 2 
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